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2D non-perturbative modeling of oscillations 
in rapidly rotating stars 
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We present and discuss results of a recently developped two dimensional non-perturbative method to compute accurate 
adiabatic oscillation modes of rapidly rotating stars . The 2D calculations fully take into account the centrifugal distorsion 
of the star while the non-perturbative method includes the full influence of the Coriolis acceleration. These characteristics 
allows us to compute oscillation modes of rapid rotators - from high order p-modes in (5Scuti stars, to low order p- and 
g-modes in /? Cephei or Be stars. 
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1 Introduction 

It is a well known fact that rotation plays a key role in 
stellar evolution. From early stages of stellar formation to 
the final steps of evolution, rotation generates various dy- 
namical processes such as meridional circulation and ro- 
tationnally induced turbulence which drive chemical ele- 
ment mixing and tra nsport of angular momentum (for ex- 
ample |Maede2|2009|). These processes are not fully under- 
stood and are still poorly modeled, but asteroseismology 
can provide important constraints provided the effects of 
rotation on stellar pulsation are better understood. In ro- 
tating stars, the centrifugal acceleration breaks the spheri- 
cal symmetry - causing distorsion - and the resonant cav- 
ity of the modes is modified. The Coriolis acceleration en- 
ters the equation of motion and affects the motion of the 
waves and the frequencies of the normal modes. For slow 
rotators, the effects of rotation on oscillation frequencies 
have be en extensively investigated with perturbative meth- 
ods (see Dziembowski & Goodcl992, Gough & Thompson 
[l990.,SaiQ.198L Soufi et al.. 1998. and references therein). 
In this approach, the angular rotation velocity Q is con- 
sidered as small compared to the oscillation frequencies, 
thereby allowing their expansion as a power series in D.. 
Perturbation methods cease to be valid whenever the rota- 
tion frequency is no longer negligible in front of the break- 
up frequency ( ^JGM/R^) or the oscillation frequency. Then 
for moderate to rapid rotators a non-perturbative treatment 
is necessary. In the non-perturbative approach, the pulsation 
equations are projected onto the spherical harmonic basis. 
The effects of the Coriolis acceleration and the stellar distor- 
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sion cause a coupling between the different spectral compo- 
nents, and the eigenvalue problem must be solved directely 
by a two-dimensionnal method. Such an approach has been 
applied to g and r-modes for uniformally rotatin g stars under 
the Cowling approximation ( ILee & Sai oI|I987|k for acoustic 
modes in unifortnally rotating polytropes ( Lignieres et al.l 
2004 iReese et al J2OO6I) in un iformally rotating ZAMS mod-| 
els jLovekin & Deupreel2008l) . and in differenti ally rotating 
ZAMS models with a co nservative rotation law ( Lovekin et al.^ 
20091 iReese et al .l2009h . We present here a 2D non-perturbative§ 
code which allows us to calculate adiabatic oscillations for 
the whole frequency range from high order g-modes to high 
order p-modes. Particular care has been taken so as to be 
able to compute pulsations for all types of stellar models. 
In the oscillation code, no hypotheses have been made on 
the fluid microphysics - non polytropic, non barotropic - 
the rotation profil is free - not necessarily conservative - 
it can be differential in radius and in latitude. The paper is 
organised as follows: in the next section, the formalism is 
explained. In section 3, we describe the numerical method 
used to solve the eigenfunction problem. A conclusion and 
perspectives follow. 



2 The formalism 



For the computation of pulsations, the stellar structure is 
reduced to its dynamical behavior. We compute oscillation 
modes as the adiabatic response of the structure to small 
perturbations - i.e. of the density, pressure, gravitation nal 
potential and velocity field - using the eulerian formalism. 
The velocity field of the equilibrium structure is only due to 



®W]LEY 

JiTiterScieTice* 



© 2010 WILEY-VCH Verlag GmbH& Co. KGaA, Weinheim 



2 



R-M. Ouazzani: 2D Non Perturbative Modeling of Oscillations 



solenoidal rotation: 



where 



vo = nxr (1) 

£1 = Q(r, 6) 008(6/) Cr - Q(r, 6) sin(0) eg (2) 



Cr and eg being the classical spherical basis vectors. Then 
the equations describing the oscillations of a self rotating 
fluid are the perturbed equations of motion: 



Idv' \ 
Po + (vo.V)v' + (v'.V)vo + p'(vo.V)vo 



= -Vp' -p'VcDo -poVO' 
The linearised continuity equation: 



d d 

^+"^|p' + V.(pov') = 



(3) 



(4) 



The Poisson equation for the perturbed gravitational poten- 
tial: 

<D' = 47TGp' (5) 
and the perturbed adiabatic relation: 

(^4)l^-fTs;)*"'(""«-TT"-'»)='' 

(6) 

By adding an auxiliary equation for the derivative of the 
gravitationnal potential t/O' - d'l>' /d(, the problem is re- 
duced to a set of first order differential and algebraic equa- 
tions. Together with specific boundary conditions, we get 
an eigenvalue problem, the eigenvalues of which are the os- 
cillation frequencies, and the eigenfunctions of which are 
the eulerian perturbations of density, pressure, gravitational 
potential, its radial derivative, and the three spatial compo- 
nents of the perturbed velocity field. 

2.1 Spheroidal geometry 

Due to the distorted shape caused by rotation, we chose to 
use a spheroidal coordinate system. This system (found by 
Bonazzola et al.lll998i) is conveni ent for setting up p roper 



boundary conditions. As done in iReese et al.l (1200a) . ( is 
defined as the radial coordinate, and is related to the spheri- 
cal r coordinate by: 

- In the stellar interior ^ 6 [0; 1], domain V: 

- 



r(f,0) = (l-e)f + 



■ {RM - 1 + e) (7) 



- In the outer domain V2, ( e [1; 2]: 

r(^,0) = 2e + (1-eK 

+ (2^^-9^2 + 12^-4)(R,(0)-l+e) (8) 

where is the colatitude, e - I - Rpoi/Req the flatness and 
Rs(0) the radius at the surface. With this mapping, the sur- 
face of the star is given by f = 1, which is very convenient 
for avoiding discontinuities at the stellar surface. At the cen- 
ter, surfaces of constant ( tend to be spherical, so that cen- 
tral regularity conditions are simplified. In the outer region, 
iso-^ surfaces regain a spherical shape at ^ - 2. 
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Fig. 1 Coordinate system used in computing the pulsation 
modes. The domain V corresponds to the star itself. V2 en- 
compasses the star, its outer limit being a sphere of radius 
r = 2 (twice the equatorial radius). 

2.2 Boundary conditions 

In order to complete the eigenvalue problem defined by the 
four equations Eq.(l3]l, Eq.©, Eq.Q, Eq.©, it is necessary 
to specify a number of boundary conditions. 
At the center of the star, where the mapping is almost spher- 
ical, the requirement is that the velocity field components 
are regular. This condition is easily expressed for functions 
expanded over the spherical harmonic basis. Concerning the 
scalar quantities, we ensure that the radial component asso- 
ciated with the spherical harmonic Y™ behaves like / (see 
Sect.3.2). 

At the surface of the star, the boundary condition is satisfied 
by means of the stressless condition 6p' - 0, which is easily 
expressed on the iso-^ surface ^ = 1 . 
It is also necessary to impose a condition on the gravita- 
tional potential so as to ensure that it goes to zero at infinity. 
With the mapping described in the former section, we can 
safely impose the classical condition on the different har- 
monic components of cD' on the outer border of V2, at ^ = 2. 

2.3 Change of variable 

In order to have solutions with a good behavior at the sur- 
face of the star, we choose to use tt' = p'/po rather than the 
pressure perturbation directly. 

3 Numerical method 

In order to isolate the radial, spheroidal and toroidal com- 
ponents of the fluid's motion, we take the radial componant 
of the equation of motion (EqO, its divergence and the ra- 
dial part of its curl. We then get a set of 7 equations and 7 
unknowns which are: the perturbed velocity vector, (n'), the 
density perturbation {p'), the perturbation of the gravitation- 
nal potential (<t)') and its derivative (d<l>7df), along with the 
related boundary conditions. 
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3.1 Projection onto the spherical harmonics 



3.2 Radial behavior 



To solve this eigenvalue problem, we develop all the vari- 
ables onto a serie of spherical harmonics, following Rieutord^ 
(fl987i) : 

For the velocity field perturbation: 



In order to ensure a proper regular behavior at the center of 
the star, and to avoid numerical convergence problems, we 
scale the radial components of the spectral decomposition 
by the appropriate powers of (. If we suppose a priori the 
regularity of the scalar variables f f (tt^, and O^), following 
Nikiforov and Uvarov (1983), these radial functions satisfy: 



(9) 



l2>\m\ 



feiO ~ C when f ^ 



(14) 



ii^ / c^Y'"(«,0) m \ 

t2>\m\ ^ ' 

^ r . m _ . 5Y-(0,<^) 

v;(^,0,<^) = - 2^ 



v', (O Y^(0, S) + w; (O- „ 

'^^^^ sine ^' 30 



where {2p - t-i + (-!)''■ For any scalar variable (tt', p', O', 



For the velocity field, the same treatment gives: 

u'^ ~ f^"', v'^ ~ if-^ and w^ ~ C,^ (15) 
This leads to the following scaling: 

'l^=^';f/ (D; = p'l^C'pt' dO^ = f'-'dC)/ 

; = ^^"'vV w^ = ^^w/ 

(16) I 



f'{C,0,ip)^ ^fliO-'^TS^A) (10) 3.3 Final eigenvalue system 



'2>|m! 

where Y" is the spherical harmonic of degree ( and az- 
imuthal order m, and f^(^) the radial functions that are to be 
determined. We include these spectral developments Eq.(|9) 
and Eg. ([Tot into the equations system Eq.(l3ll to Eq.®. We 
get a linear partial differential equations system in terms of 
the variables u^ , v^^, wj,^, n'^^, c/<tj,^, and p^^ which can 
formally be written as: 



^ e{y^I{0,4>)) = 

f2>|m| 



(11) 



In the general case - where rotation breaks the spherical 
symmetry - the problem is not separable in (, and 0. As a 
result, a coupled equations system is obtained by projecting 
Eq.lfTTI) onto the spherical harmonics basis: 



> Im|, ^ 

^2>|m| 



+ O0 ^ 



sinfded^ 



E(Y^)Yr 



(12) 



In the set of equations, 4 are differential equations for the 
variables m', tt', O' and c/ft', with respect to the radial coor- 
dinate ^, and 3 are not: 

^ = (All + 5o-Ai2)yi + (A21 + 60- A22)y 2 (17) 

= (Bii + 5o-Bi2)yi + (B21 + 5o-B22)y2 (18) 
where cr — ctq + bcr. y\ and y2 are the column vectors - 
with 4 M and 3 M components respectively - containing 
the unknown coefficients of the spectral decomposition: 

yi = (7f^^,...,7r;„,dO;,,...,dO;„,<l)^_,..., 0)^^,0^^, ...,u^^) 

y2 = (v;,,...v^„,w^_,...w;^,p;^,...p^j 

where 

^1 = |m|+p and = |m| + 2(M - 1) + p (19) 

Using ( fTSl ). a matrix inversion allows to express y2 as a 
linear function of yi . The system is then reduced to a differ- 
ential system of 4 independent variables: 



where Y™* is the complex conjugate of Y™. Finally, to get 
a finite and equal number of equations and variables, we 
truncate the series at the 2M''' term. Given that the equilib- 
rium model is symmetrical with respect to the equator, half 
of the terms are dropped in the projection, and the system 
only couples terms with the same symmetry, i.e. parity. The 
selection rules operating here are: 

^ = |m|+2(j-l) + p, j6[l:M] (13) 

p = if m and t are of the same parity, 
p = 1 otherwise. 

We then are able to solve the eigenvalue problem separately 
for the even and odd eigenfunctions, and for a given az- 
imuthal order 



^ -(A + 5tr.A,Jyi 



3.4 Radial resolution: finite differences 



(20) 



Considering two consecutive layers i and i+1, such that 
^(i + 1) - 4"(i) = h, a Taylo r development of any function y 
gives dScuflaire et al.ll2008b : 



h dy 



h2 d^y 



2 d?« ^ 12 d?® = 



12 d^2 



We can develop a 5* order finite differences scheme that 
only involves two consecutive layers (f, and f,+i). 
We then obtain a global eigenvalue system: 



AAY ^^crAA^^Y 



(22) 
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Where AA and AA^o- are block diagonal matrices. The block§ 
i couples the layers i and i+1 and: 



Y = 



yi(N). 



(layers 1, 



N) 



(23) 



3.5 Inverse iteration algorithm 



In order to solve Eg. (122b w e use a genera hzation of the in- 
verse iteration method (see lDupretll200lb . Starting from a 
first estimate of the eigenvector Yq and eigenfrequency cor- 
rection 6o-Q (we take dcro - 0, i.e. cr = ctq), we compute the 
next step in the iteration using the formula: 



Yk+i = AA-'AA*,Yk 



(24) 



provided AA"' AA^o- is diagonalizable. Solving the eigen- 
value problem is then equivalent to solving the linear sys- 
tem: 



AAYk+i = AAa^Yk 



(25) 



To do so, we perform a LU factorization of the matrix AA - 
L U, where L is a lower triangular matrix and U an upper 
triangular one. Therefore, the LU factorization needs to be 
made only once, and the only thing to be done afterwards is 
to solve at each step of the inverse iteration the two triangu- 
lar systems: 



LX = AAs^Yk 
UYk+i = X 



(26) 



To avoid ill-conditionning problems, we adopt a special kind§ 
of pivoting strategy where the pivoting is done alternatively 
on the columns and on the lines of the matrix. Thanks to 
the chosen finite difference resolution (Sect. 3.4), AA and 
AA^o- are block diagonal matrices, therefore, L and U are 
also block diagonal matrices, where the blocks are triangu- 
lar. The non-zero elements are all located inside the blocks 
and, during the algorithm, the permutations between lines 
and the permutations between columns will keep the non- 
zero elements in the same block. This allows us to keep the 
narrow band shape for the system, and reduces the memory 
needed and the computational time. Once the eigenvector Y 
is computed with sufficiently high precision, the eigenvalue 
can be calculated using a generalization of the Rayleigh ra- 
tio: 



6cr 



Y* AA'_^AAY 
Y* AA*^AA5^Y 



(27) 



where Y* and AA^^ are the hermitian conjugate of Y and 
AA^o-. The value of 6cr given by Eq. dZTl i minimizes: 

S- ^\ (AA - 6(TAAsa)Y \\- (28) 

The solution is then obtained when a certain iteration cri- 
terium is reached, which depends on the precision we want 
on the oscillation frequency. 



4 Conclusion and perspectives 

We presented here a new code that performs 2D non-perturbative§ 

calculations of adiabatic oscillations for all kinds of stellar 
structures - stratified, differentially rotating. We used a nu- 
merical method which efficiently saves computationnal time 
and memory, and would make the code available for seismic 
interpretation. It is now under a series of test - comparison 
with perturbative methods for evolved stellar models, and 
comparison with non-perturbative computations for poly- 
tropic models. With this new tools, we aim at modeling os- 
cillations of 2D stratified models of stars. A first step will be 
to take as the equilibrium model a ID evolved stellar model 
where the whole effect of rotation on the microphysics have 
been ta ken into account, and use a self-consistent method 
(used in lRoxburghll2006[ for example) in order to compute 
the distorsion due centrifugal acceleration. 
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